{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e16fe530-d09a-451c-a888-c9b7d177c093",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "ReadOnlyMemoryError()",
     "output_type": "error",
     "traceback": [
      "ReadOnlyMemoryError()",
      "",
      "Stacktrace:",
      " [1] top-level scope",
      "   @ In[9]:197"
     ]
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "using Random\n",
    "using Plots\n",
    "using Statistics \n",
    "using DataFrames\n",
    "using DelimitedFiles\n",
    "\n",
    "# 生成初始状态\n",
    "function generate(n)\n",
    "    return zeros(Int, n*n)\n",
    "end\n",
    "\n",
    "# 计算蜂窝晶格的近邻\n",
    "function neighbor_honeycomb(n)\n",
    "    nbor1 = zeros(Int, 3, n * n)\n",
    "    nbor2 = zeros(Int, 6, n * n)\n",
    "    nbor3 = zeros(Int, 3, n * n)\n",
    "\n",
    "    for ispin in 1:(n * n)\n",
    "        iy = div(ispin - 1, n)\n",
    "        ix = ispin - 1 - iy * n\n",
    "\n",
    "        # 计算周期性边界条件下的相邻自旋的索引\n",
    "        ixp = mod(ix + 1, n)\n",
    "        iyp = mod(iy + 1, n)\n",
    "        ixm = mod(ix - 1 + n, n)\n",
    "        iym = mod(iy - 1 + n, n)\n",
    "\n",
    "        # 计算二阶近邻的索引\n",
    "        ixpp = mod(ix + 2, n)\n",
    "        ixmm = mod(ix - 2 + n, n)\n",
    "        iypp = mod(iy + 2, n)\n",
    "        iymm = mod(iy - 2 + n, n)\n",
    "\n",
    "        # 一阶近邻\n",
    "        nbor1[1, ispin] = iy * n + ixp + 1\n",
    "        nbor1[3, ispin] = iy * n + ixm + 1\n",
    "\n",
    "        if iy % 2 == 1\n",
    "            if ix % 2 == 1\n",
    "                nbor1[2, ispin] = iyp * n + ix + 1\n",
    "\n",
    "                # 二阶近邻\n",
    "                nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                # 三阶近邻\n",
    "                nbor3[1, ispin] = mod(iyp * n + ixpp, n*n) + 1\n",
    "                nbor3[2, ispin] = mod(iyp * n + ixmm, n*n) + 1\n",
    "                nbor3[3, ispin] = mod(iym * n + ix,   n*n) + 1\n",
    "\n",
    "            else\n",
    "                nbor1[2, ispin] = iym * n + ix + 1\n",
    "\n",
    "                # 二阶近邻\n",
    "                nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                # 三阶近邻\n",
    "                nbor3[1, ispin] = mod(iyp * n + ix, n*n) + 1\n",
    "                nbor3[2, ispin] = mod(iym * n + ixmm, n*n) + 1\n",
    "                nbor3[3, ispin] = mod(iym * n + ixpp, n*n) + 1\n",
    "\n",
    "            end\n",
    "        else\n",
    "            if ix % 2 == 1\n",
    "                nbor1[2, ispin] = iym * n + ix + 1\n",
    "\n",
    "                # 二阶近邻\n",
    "                nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                # 三阶近邻\n",
    "                nbor3[1, ispin] = mod(iyp * n + ix, n*n) + 1\n",
    "                nbor3[2, ispin] = mod(iym * n + ixmm, n*n) + 1\n",
    "                nbor3[3, ispin] = mod(iym * n + ixpp, n*n) + 1\n",
    "\n",
    "            else\n",
    "                nbor1[2, ispin] = iyp * n + ix + 1\n",
    "\n",
    "                # 二阶近邻\n",
    "                nbor2[1, ispin] = iy  * n + ixpp + 1\n",
    "                nbor2[2, ispin] = iyp * n + ixp  + 1\n",
    "                nbor2[3, ispin] = iyp * n + ixm  + 1\n",
    "                nbor2[4, ispin] = iy  * n + ixmm + 1\n",
    "                nbor2[5, ispin] = iym * n + ixm  + 1\n",
    "                nbor2[6, ispin] = iym * n + ixp  + 1\n",
    "\n",
    "                # 三阶近邻\n",
    "                nbor3[1, ispin] = mod(iyp * n + ixpp, n*n) + 1\n",
    "                nbor3[2, ispin] = mod(iyp * n + ixmm, n*n) + 1\n",
    "                nbor3[3, ispin] = mod(iym * n + ix, n*n) + 1\n",
    "\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return nbor1, nbor2, nbor3\n",
    "end\n",
    "\n",
    "# 计算能量和磁矩\n",
    "function calculate(ising, nbor1, nbor2, nbor3, J1, J2, J3, N, U, μ)\n",
    "    e = 0.0\n",
    "    neighbor_sum = 0.0\n",
    "    for i in 1:N\n",
    "        for j in 1:3\n",
    "            neighbor_sum = neighbor_sum - J1 * ising[i] * ising[nbor1[j, i]]\n",
    "        end\n",
    "        for j in 1:6\n",
    "            neighbor_sum = neighbor_sum - J2 * ising[i] * ising[nbor2[j, i]]\n",
    "        end\n",
    "        for j in 1:3\n",
    "            neighbor_sum = neighbor_sum - J3 * ising[i] * ising[nbor3[j, i]]\n",
    "        end\n",
    "    end\n",
    "    e = e + neighbor_sum\n",
    "    e = e / 2\n",
    "\n",
    "    m = sum(ising)\n",
    "    for i in 1:N\n",
    "        e = e + (U / 2) * ising[i] * (ising[i] - 1)\n",
    "    end\n",
    "\n",
    "    for i in 1:N\n",
    "        e = e - μ * ising[i]\n",
    "    end\n",
    "\n",
    "    return e, m\n",
    "end\n",
    "\n",
    "# 生成坐标\n",
    "function generate_coordinates(N, M)\n",
    "    x_coords = zeros(N, M)\n",
    "    y_coords = zeros(N, M)\n",
    "\n",
    "    for i in 1:N\n",
    "        if i % 2 == 1\n",
    "            base_value = (i - 1) * 1.5\n",
    "        else\n",
    "            base_value = (i - 2) * 1.5 + 1\n",
    "        end\n",
    "\n",
    "        for j in 1:M\n",
    "            x_coords[i, j] = (j - 1) * sqrt(3) / 2\n",
    "\n",
    "            if i % 2 == 1\n",
    "                y_coords[i, j] = j % 2 == 1 ? base_value : base_value - 0.5\n",
    "            else\n",
    "                y_coords[i, j] = j % 2 == 1 ? base_value : base_value + 0.5\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    return x_coords, y_coords\n",
    "end\n",
    "\n",
    "n = 4\n",
    "N = n * n\n",
    "\n",
    "# 定义参数范围\n",
    "J1_values = [1]\n",
    "J2_values = [-0.5]\n",
    "J3_values = [0.1]\n",
    "U_values = 0.7:0.2:0.9\n",
    "μ_values = -0.1:0.2:0.3\n",
    "\n",
    "for J1 in J1_values\n",
    "    for J2 in J2_values\n",
    "        for J3 in J3_values\n",
    "            for U in U_values\n",
    "                for μ in μ_values\n",
    "                    # 创建保存图形的文件夹\n",
    "                    save_dir = \"results_J1_$(round(J1, digits=2))_J2_$(round(J2, digits=2))_J3_$(round(J3, digits=2))_U_$(round(U, digits=2))_μ_$(round(μ, digits=2))\"\n",
    "                    ispath(save_dir) || mkpath(save_dir)\n",
    "\n",
    "                    e_steps = zeros(3^N)\n",
    "                    m_steps = zeros(3^N)\n",
    "                    ising_steps = zeros(3^N, N)\n",
    "                    ising = generate(n)\n",
    "                    nbor1, nbor2, nbor3 = neighbor_honeycomb(n)\n",
    "\n",
    "                    for k1 in 0:3^N - 1\n",
    "                        num1 = lpad(string(k1, base=3), N, '0')\n",
    "                        for j in 1:N\n",
    "                            ising[j] = parse(Int, num1[j])\n",
    "                        end\n",
    "\n",
    "                        ising_steps[k1+1,:] = ising[:]\n",
    "                        e_steps[k1 + 1], m_steps[k1 + 1] = calculate(ising, nbor1, nbor2, nbor3, J1, J2, J3, N, U, μ)\n",
    "                    end\n",
    "\n",
    "                    min_energy, index = findmin(e_steps)\n",
    "                    count_min_energy = count(x -> x == min_energy, e_steps)\n",
    "                    indices = findall(x -> x == min_energy, e_steps)\n",
    "                    min_configurations = ising_steps[indices, :]\n",
    "\n",
    "                    x_coords, y_coords = generate_coordinates(n, n)\n",
    "                    x = reshape(x_coords', :)\n",
    "                    y = reshape(y_coords', :)\n",
    "\n",
    "                    for i in 1:count_min_energy\n",
    "                        p = plot()\n",
    "                        z = min_configurations[i,:]\n",
    "                        # 修改颜色映射规则\n",
    "                        colors = [\n",
    "                            if z_i == 0\n",
    "                                :black\n",
    "                            elseif z_i == 1\n",
    "                                :blue\n",
    "                            elseif z_i == 2\n",
    "                                :red\n",
    "                            else\n",
    "                                :black # 处理其他可能的值，默认设为黑色\n",
    "                            end\n",
    "                            for z_i in z\n",
    "                        ]\n",
    "                        scatter!(x, y, color=colors, aspect_ratio=1)\n",
    "\n",
    "                        # 保存图形\n",
    "                        savefig(p, joinpath(save_dir, \"configuration_$(i).png\"))\n",
    "                    end\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef7a903-cd93-4845-a8f7-b54d04669dbc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.4",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
